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We investigate the behavior of the slip length in Newtonian liquids subject to planar shear bounded 
by substrates with mixed boundary conditions. The upper wall, consisting of a homogenous surface 
of finite or vanishing slip, moves at a constant speed parallel to a lower stationary wall, whose 
surface is patterned with an array of stripes representing alternating regions of no-shear and finite 
or no-slip. Velocity fields and effective slip lengths are computed both from molecular dynamics 
(MD) simulations and solution of the Stokes equation for flow configurations either parallel or 
perpendicular to the stripes. Excellent agreement between the hydrodynamic and MD results is 
obtained when the normalized width of the slip regions, a/o > O(10), where a is the (fluid) molecular 
diameter characterizing the Lennard- Jones interaction. In this regime, the effective slip length 
increases monotonically with a/a to a saturation value. For a/o < O(10) and transverse flow 
configurations, the non-uniform interaction potential at the lower wall constitutes a rough surface 
whose molecular scale corrugations strongly reduce the effective slip length below the hydrodynamic 
results. The translational symmetry for longitudinal flow eliminates the influence of molecular scale 
roughness; however, the reduced molecular ordering above the wetting regions of finite slip for small 
values of a/o increases the value of the effective slip length far above the hydrodynamic predictions. 
The strong correlation between the effective slip length and the liquid structure factor representative 
of the first fluid layer near the patterned wall illustrates the influence of molecular ordering effects 
on slip in non-inertial flows. 

PACS numbers: 61.20.Ja, 68.08.-p, 68.35.Af, 83.10.Rs, 83.50.Rp, 83.60.Yz 



I. INTRODUCTION 

The development of micro- and nanofluidic devices for 
the manipulation of films, drops and bubbles requires 
detailed knowledge of intcrfacial phenomena and small 
scale flows. These systems, which are distinguished by a 
large surface-to- volume ratio and flow at small Reynolds, 
capillary and Bond numbers, are strongly influenced by 
boundary effects p). Liquid affinity to nearby solid 
boundaries can be reduced through chemical treatments 
|H [1 LULl substrate topology |a0> 13 or the nucleation 
of nanobubbles on hydrophobic glass surfaces [jj 0, 0] . 
Weak van der Waals interactions between a polymer melt 
and solid wall 0, 0, 0] or between two immiscible 
polymers can also lead to significant slippage and 
reduced frictional resistance. The degree of slip is nor- 
mally quantified through the slip length defined as the 
distance from the surface within the solid phase where 
the extrapolated flow velocity vanishes [ly. Numerous 
experimental and theoretical studies have examined how 
the slip length is influenced by such factors as the de- 
gree of hydrophobicity la, [l7l. the substrate to p ogr aphy 
and surface roughness [& Hill 111 IM EJ lM pi pi , 
the presence of interstitial lubricating l ayer s [221 Y2A [25j , 
the polymer molecular we ight HiL 1261 |27| and the ap- 
plied shear rate 0, |H |H [3(], l3lj . In a recent develop- 
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ment, the large values of the slip length extracted from 
experiments involving the pressure-driven flow of water 
through hydrophobically coated capillaries have been at- 
tributed 32, 33J to the spontaneous nucleation of a dense 
and stable layer of nanobubbles in water films adjacent 
to hydrophobic glass surfaces jjj El 0] . Of special in- 
terest is the corresponding reduction in drag achieved 
by proportional substitution of liquid-solid contact area 
with liquid-gas contact area or equivalently, substitution 
of regions of no-slip or finite slip by regions of essentially 
no-shear (i.e. infinite slip). 

Interest in the hydrodynamic behavior of liquid films 
in the vicinity of surfaces with mixed boundary condi- 
tions dates back several decades to the work of Philip 
[34l |3f| . He examined the steady flow of an incompress- 
ible and inertia-less Newtonian liquid driven either by 
a uniform shear stress or uniform pressure gradient and 
subject to mixed wall boundary conditions. These were 
represented by surfaces consisting of alternating striped 
regions of no-shear and no-slip among other geometries. 
Using conformal mapping, Philip [34 . l35| derived ana- 
lytic expressions for the streamfunction and volumetric 
flux for flow perpendicular (transverse configuration) or 
parallel (longitudinal configuration) to the striped ar- 
ray in the limit of Stokes flow. Recently, Lauga and 
Stone [3j| investigated the behavior of the effective slip 
length for steady Poiseuille flow through a capillary of 
circular cross-section whose inner wall consists of period- 
ically distributed regions of no-slip and no-shear. Philip's 
earlier treatment was used to extract the slip length for 
longitudinal configurations; additional analysis was re- 
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FIG. 1: (Color online) (a) Transverse and (b) longitudinal 
flow orientations for a liquid film subject to planar shear in a 
cell with wall separation d. Darker stripes of width a signify 
regions of finite slip or no-slip. White stripes signify regions of 
no-shear (or equivalently perfect slip). The upper wall moves 
at constant speed U relative to the lower stationary surface 
(z = 0). The periodicity of the lower wall pattern geometry is 
designed by A. 

quired for transverse configurations. Comparison of their 
results with available experimental measurements sug- 
gests what model parameter values would reproduce the 
experimental slip lengths. For slip lengths in the nanome- 
ter range, one might ask whether a hydrodynamic anal- 
ysis can correctly predict these values or whether the 
molecular aspects of the fluid can strong influence the slip 
behavior causing deviations from the continuum theory. 

Molecular dynamics (MD) simulations provide an ideal 
tool for investigating the conformation and behavior of 
fluid molecules adjacent to chemically or topologically 
textured substrates. The boundary conditions which es- 
tablish the flow profile are not specified a priori, but arise 
naturally from the wall-fluid contrast in density and the 
fluid-fluid and wall-fluid interaction potentials. In recent 
years, many groups have examined how various molec- 
ular parameters characterizing the wall and fluid prop- 
erties affect the degree of slip at liquid-solid interfaces. 
In particular, it has been demonstrated that the struc- 
ture factor and contact density representative of the first 
fluid layer adjacent to a wall significantly influence the 
deg ree of slip in Newtonian and non-Newtonian fluids 
H HI IM |M HI H H3. The results of this current 
study confirm the importance of these molecular param- 
eters for flow on heterogeneous substrates. 

In this work, we investigate the behavior of the slip 
length in viscous films under planar shear bounded by 
substrates with mixed boundary conditions using both 
molecular dynamics (MD) simulations and Stokes flow 
computations. The upper wall, consisting of a homoge- 
nous surface of finite or no-slip, moves at a constant 
speed, U, a distance d above a lower stationary wall, 
whose surface is patterned with an infinite array of stripes 
representing alternating regions of no-shear and finite 
or no-slip. As shown in Fig. ^ we consider transverse 
and longitudinal flow configurations and compute the 
corresponding velocity fields and effective slip lengths 
for a wide range of stripe widths, periods and liquid- 
solid affinities. Excellent agreement between the hy- 
drodynamic and MD results is obtained when the nor- 
malized width of the slip regions, a/a > O(10), where 



a is the (fluid) molecular diameter characterizing the 
Lennard-Jones interaction. For surface patterns ap- 
proaching molecular size, the degree of fluid ordering near 
the patterned wall, as quantified by the in-plane struc- 
ture factor and contact density in the first liquid layer, 
plays a dominant role causing significant deviations from 
the hydrodynamic predictions. These deviations can be 
explained in the context of effective surface roughness 
and molecular ordering effects. 

II. HYDRODYNAMIC ANALYSIS 

In the limit of vanishingly small Reynolds number 
Re = pUdj '(j,, where p and p denote the (constant) liquid 
density and viscosity, inertial effects are negligible. The 
velocity profile is then governed by the Stokes equation, 
V 2 u = Vp//i, where the velocity field, u, satisfies the 
condition of incompressibility, Vu = 0, andp denotes the 
pressure distribution which in this study is induced by the 
patterned substrates. Application of the divergence oper- 
ator to the Stokes equation shows that the pressure field 
satisfies the equation V 2 p = 0. It then follows that the ve- 
locity field satisfies the biharmonic equation V 2 V 2 u = 

m 

In the next section, we derive the boundary conditions 
(BCs) corresponding to transverse [Fig. Qfa)] and longi- 
tudinal [Fig.njb)] flow orientations. These conditions are 
used to compute numerical solutions of the streamfunc- 
tion, velocity field and effective slip length as a function 
of the dimensionless stripe width of the finite slip regions, 
a/A, and the dimensionless surface period X/d. The y- 
axis is oriented parallel to the stripe edges for cither 
configuration. All numerical calculations were performed 
with the finite element software FemLab 2.3 |4l],[42] using 
triangular elements with quadratic basis functions. The 
solutions reported converged upon mesh refinement. 



A. Transverse configuration 

The two dimensional velocity field corresponding to 
the transverse configuration shown in Fig. ^a) is rep- 
resented by u(x,z) = (u,0,w) — (dip/dz, 0, — dip/dx), 
where ip(x, z) denotes the streamfunction, which implic- 
itly satisfies the continuity equation V-u = 0. The vortic- 
ity vector f2 = Vx u= (0, u>, 0), where u> — du/dz—dw/dx, 
has only one non-zero component. According to these 
definitions, it follows that 

1. Boundary Conditions 

Solutions of the equations for the vorticity and stream- 
function given by Eq. require the specification of eight 
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BCs. The computational domain sketched in Fig.JSfa) is 
defined by the region bounded by the upper and lower 
walls (0<2:<d) and the dashed lines (0<x<A/2) corre- 
sponding to the midplanes of neighboring stripes. White 
surfaces designate shear-free boundaries (i.e. surfaces of 
perfect slip) ; dark surfaces designate boundaries of finite 
or no-slip. Throughout, partial derivatives are denoted 
by letter subscripts e.g. dtjj/dx = t\j x . 

The top and bottom walls represent impenetrable sur- 
faces where w(x, z = 0) = w(x 1 z = d) = 0, or in terms of 
the streamfunction, ip x (x, z = 0) = ip x (x, z = d) = 0. The 
tangential component of the velocity field must satisfy 
mixed slip and shear conditions at the upper and lower 
walls of the cell. The no-shear BC is given by u z (0 < 
x <- A^a ^ z = o) — 0. Slip surfaces are characterized by 
the Navier [l(| slip condition u <x< | , z = 6\ = bu z 
and u[x, z = d] = U — bu z . The Navier slip length b is as- 
sumed constant, i.e. independent of the shear rate 7. 

The lateral boundary conditions for the scalar field 
u are derived from the following symmetry considera- 
tions. The biharmonic equation V 2 V 2 u = involves 
x-derivatives of even order only. The lower wall com- 
prises an infinite number of mirror symmetry planes lo- 
cated at the at the stripe centers x = nX/2, for all in- 
tegers n. Since the upper surface is homogeneous and 
translationally invariant, the mirror symmetry imposed 
by the lower surface determines which symmetry applies 
throughout the entire Couette cell. The scalar field, u, 
therefore also assumes mirror symmetry about the stripe 
centers such that u(x, z) = u(—x, z) and u x (x, z) = for 
all x — nX/2 and integers n. From the continuity equa- 
tion, it then also follows that w z (x = n\/2) = 0, i.e. w is 
independent of the coordinate z within any mirror plane. 
Since the upper and lowers walls are impenetrable, i.e. 
w(x, z = 0) = w(x, z = d) — 0, this constraint reduces to 
the BC 4> x (x = 0) = = ip x (x = \/2). 

The continuity equation requires u x + w z — 0. To- 
gether with the condition u(x,z) = u(—x,z), this im- 
plies w(x, z) = —w{—x, z) such that w(x, z) = and 
w xx {x, z) = at all x = n\/2 where n = 0, 1, 2 . . .. Sub- 
stitution of this last relation and u x (x = 0, z) — into 
the expression for u> x leads to to x (x — nX/2, z) — 0. Re- 
gions of no-shear at the lower wall are represented by the 
condition lu(x, z — 0) = 0. 

Along the top and bottom walls, the scalar compo- 
nent w is independent of the coordinate x and there- 
fore w x (z = 0) = w x (z — d) — 0. Consequently, 
the vorticity at the top and bottom walls reduces to 
lu = u z and the Navier slip conditions can be rewritten as 
V^z (^<x<^,z = 6) = buj and ip z (x,z = d) = U - buj. 
The relation w(x, z = 0) = w(x, z = d) — ipx — also 
implies that the streamfunction is constant in the planes 
z = and z = d, whose values we denote by V'top and 
V'bottom- The difference in the streamfunction value be- 
tween the top and bottom walls is equal to the volumetric 
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FIG. 2: (Color online) Computational domain and bound- 
ary conditions used for solution of the Stokes equation corre- 
sponding to the (a) transverse and (b) longitudinal flow ori- 
entation shown in Fig. The computational domain consists 
of the region bounded by the upper and lower walls (z = and 
z = d) and the lateral dashed lines (a; = and x = X/2), which 
are positioned at neighboring midplanes of the no-shear and 
finite-slip regions. 



flux per unit length along the j/-axis [40| since 
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Because the streamfunction can only be determined 
within an arbitrary constant, we set the value of ^bottom 
in these studies to zero without loss in generality. The 
complete set of BCs for the vorticity and streamfunctions 
are therefore given by: 
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(3) 
(4) 
(5) 
(6) 
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2. Solution procedure 



The value of tptop is determined from the pressure field 
as follows. The Stokes equation for the vertical compo- 
nent of the velocity field is given by w xx + w zz =p z / fi. As 
argued in a previous section, however, w xx (x = nX/2, z) = 
0, and since w is independent of z along any mirror sym- 
metry plane, w zz (x — nX/2, z) — 0. The pressure is there- 
fore independent of the vertical coordinate z in all planes 
x = nX/2. Furthermore, in the absence of any externally 
applied pressure gradient, as is the case here, and be- 
cause of the flow periodicity, p(x = 0) = p(x = A). Since 
it was previously argued that u exhibits mirror symme- 
try about the planes x — nA/2, it must also be true of 
p x since u xx + u zz — p x j [i. Consequently, the pressure 
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is equal at the lateral boundaries of the computational 
cell, i.e. p(x = 0) = p(x = A/2). For convenience we set 
p(x — 0) = 0. This constraint, coupled with the relation 
Px/ ^ = — (V x SI) -e x = u> z , was used to adjust the numer- 
ical value of iptop by requiring that the following integral 
vanish identically: 



A/2 p\/2 

p x dx — /i I u z dx = 

I! JO 



(10) 



The value of the effective slip length, L s , corresponding 
to the overall flow within a patterned cell was obtained 
from linear extrapolation of the averaged velocity profile 

A/2 

(u) = (2/A) J u(x,z)dx to zero. Since at planes of 
mirror symmetry, u x = and p(x = 0) = p(x = A/2), the 

integral n J^ 2 V 2 u(x, z)dx — j^ 2 p x dx — 0, reduces to 
(u)zz — 0. The averaged velocity field, (u), is therefore a 
linear function of z and geometric similarity establishes 
the relation for the effective slip length, namely 
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(u(z = 0)) 
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(11) 



For the numerical analysis, the equations for the vorticity 
and streamfunction given by Eq. and the BCs given 
by Eqs. (|3I9[1 were non-dimcnsionalized according to the 
rescaled variables 
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leading to 
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In Section III A 41 we present numerical solutions to 
Eqs. (|15fl and the extracted values of L s as a function 
of the local slip length b and pattern geometry. Ana- 
lytic expressions are derived in the limits X/d — > and 
X/d — > oo. 



3. Perturbative analysis for b — 

In order to enhance the numerical precision of solu- 
tions corresponding to small values of L s , the velocity 
and pressure fields were decomposed into two contribu- 
tions, u = u + ui andp=po+Pi- Here, Uo = (Uz/d, 0, 0) 
and po — correspond to the velocity and pressure fields 
for planar shear flow subject to no-slip at both solid 
boundaries. The Stokes equation then reduces to the 
form /xV 2 Ui = Vpi, where the perturbed velocity field 
satisfies the continuity equation V-Ui =0. The follow- 
ing BCs for the perturbed streamfunction and vorticity 



fields were determined in similar fashion as those in Sec- 
tion EXU 
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where tp±^ =U\, ipx, x = ~wi and u)\ =iti, z — w±, x . Non- 
dimensionalization of the vorticity and streamfunction 
perturbations toi and ipi as in Section Hi A 21 leads to: 
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4- Numerical results and limiting cases 



In Fig. |3{a) is plotted the numerical results for the 
normalized effective slip length, L s /d, as a function of 
the aspect ratio, X/d, for the transverse configuration. 
Over the range shown, L s /d increases monotonically with 
X/d, saturating at a constant value beyond X/d ~ O(10). 
When X/d — > oo, any significant variation in the veloc- 
ity and pressure fields will be localized near the plane 
x = (A — a)/2, where the BCs change from no-shear 
to finite slip. Since p(x = 0) = p(x = A/2), the longi- 
tudinal average of the lateral pressure gradient within 
the cell must vanish (i.e. (p x ) = 0) and any pressure 
gradient above the surface of no-shear will be canceled 
by an opposing gradient above the surface of finite slip. 
Since the transition region in the vicinity of the the plane 
x = (A— a)/2 does not contribute significantly in the limit 
X/d^oo, the condition (p x } = is equivalent to the con- 
dition 
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where the subscripts 1 and 2 refer to the regions above 
the surface of no-shear (1) and finite slip (2). 

Now we first consider the case b/d — 0. Since the flux 
must remain constant, 



u\dz = I uidz 
Jo 
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It follows that (px)2 = &{Px)i — Q/iU/d 2 , which when 
coupled with Eqs. and (|25[) . yields the limiting value 
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The same analysis can be extended to the case b/dy^O 
with the general result: 



lim — p 

X/d— too d 



Xd 2 + 8Xbd + 12Xb 2 - ad 2 - Aabd 
Aad{d + 36) 



(30) 



The horizontal asymptotes (dotted lines) shown in 
Fig. |^a) for X/d > 10 represent solutions to Eq. l|3Tj|) 
for the designated values of b/d and a/X. 

In the opposite limit X/d — > 0, i.e. where the upper 
and lower walls are essentially infinitely far apart, the 
deviation of the flow field from pure shear flow over a 
homogeneous surface with slip length L s is limited to a 
thin layer whose thickness scales with A. As a conse- 
quence, the effective slip length should be independent 
of the cell depth d and independent of the particular 
mechanism used to generate the flow, i.e. the same slip 
length should result for pressure-driven or shear-driven 
flow. Lauga and Stone 33] determined the asymptotic 
behavior of the effective slip length for pressure-driven 
flow in a cylindrical tube of radius R with periodically 
distributed (transverse) rings denoting alternating sur- 
faces of no-shear or no-slip (6 = 0): 
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The solutions to Eq. H31|) , obtained by replacing the cap- 
illary radius R with the cell depth d, superimpose per- 
fectly (sloped dashed lines) onto the full numerical solu- 
tions shown in Fig. Ol^a). In this limit, the slip length 
increases linearly with X/d up to a limit A/d ~ 1. 

The effective slip length in the limit X/d — > for the 
case 6^0 can be derived as follows. When the array 
period A is much smaller than the local slip length 6, 
the slip velocity u(x, z — 0) should saturate towards a 
constant value, u s o, over the entire interval < x < A/2. 
Since u(z = 0) = bu z (z = 0), it is also expected that the 
velocity gradient, u z (z — 0), will assume a constant value 
in the region [A — a)/ 2 < x < A/2, z — 0]. Consequently, 
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i.e. the effective slip length becomes independent of X/d 
for fixed a/A. The term proportional to (A— a) /A accounts 
for the vanishing contribution of the no-shear regions to 
d(u)/dz(z = 0). The horizontal dashed lines shown in 
Fig. OUa) for A/d<0.1 represent solutions to Eq. lj3~2l for 
the designated values of b/d and a/X. 

In Fig. |^b) is plotted the effective slip length, L s /d, 
versus a/X for b/d — and X/d — 1.0 and 14.3. The 
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FIG. 3: (Color online) (a) Normalized slip length, L s /d, ver- 
sus normalized pattern period, X/d, derived from the Stokes 
solutions for the transverse flow orientation. The parameters 
values shown are b/d — 0, 0.048 and 0.098 and a/A = 0.25 and 
0.50. The straight lines superimposed on the numerical solu- 
tions for b/d = correspond to the analytic limit L s /d ~ X/d. 
(b) Normalized slip length L 3 /d versus normalized stripe 
width a/X for X/d = 1 and 14.3 in the limit b/d = 0. The 
dashed lines correspond to the function L 3 /d = A\n[a/X] + B, 
with fitting parameters A and B. The data points for the 
case X/d = 14.3 (A = -2.245, B = -3.952) are scaled by a 
factor 0.3 for convenience, (c-d) Streamlines corresponding 
to the transverse Stokes flow solutions for (c) A/d=l and (d) 
A/d = 20 where b/d = 0.048 and a/X = 0.5. The cell domain 
size is A/2 x d; the thin vertical lines designate the boundary 
between surfaces of no-shear (left) and finite slip (right). 



data points for X/d = 14.3 are scaled by a factor 0.3 for 
convenience. The effective slip vanishes as a/A— >1 since 
the surface coverage by regions of perfect slip decreases 
to zero. The numerical results were compared to a Taylor 
expansion of Eq. I|31|) in the limit of a/X — >0 
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The dashed line shown in Fig. &b) for X/d — 1.0 per- 
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fectly superimposes on the results of the full numerical 
solutions. The numerical solution for X/d = 14.3 can 
also be approximated by a fit-function A log [a/ A] + B for 
a/A <c 1, as shown by the dashed line; however, Eq. I)33|) 
no longer holds because X/d ^ 1. 

Streamlines of the flow field, corresponding to the con- 
tour lines (i.e. constant values) of the streamfunction, 
are shown in Figs. EIc,d). The left and right panels rep- 
resent the solutions for b/d — 0.048 and a/A = 0.5 for 
(c) X/d = 1 and (d) 20. The vertical line denotes the 
transition in boundary condition at the lower wall from 
no-shear (left) to finite slip (right). For small values of 
X/d < 1.0, the streamlines are essentially horizontal in 
the larger portion of the cell and the deviation of the 
streamfunction from pure Couette flow over a homoge- 
neous surface is confined to a small distance from the 
patterned wall. As X/d increases, the perturbation ex- 
tends further away from the lower boundary. For large 
X/d the streamlines are horizontal above the individual 
stripes except for a step-like vertical displacement in the 
vicinity of the transition point x — (A — a)/2. 



B. Longitudinal configuration 

The velocity field corresponding to the longitudinal 
configuration shown in Fig. ^b) is unidirectional and 
given by u(x,z) = (0,u,0). There is no pressure gra- 
dient in this configuration and the numerical solutions 
are derived directly from the Stokes equation V 2 u = 0. 
The computational cell is shown in Fig.^b), where the 
direction of motion of the upper wall is indicated by the 
white concentric circles. Only four BCs are required for 
solution of the velocity field v. Aside from the obvious 
constraints of finite slip, v x must vanish at x = and 
x = A/2 because these are planes of mirror symmetry. 
The complete set of BCs is given by: 

v(0<x<±,z = d) = U-bv z {x,z = d) (34) 

v z (0<x<^,z = 0) = (35) 

v(^<x<±,z = 0) = bv z (x,z = 0) (36) 

v x (x=0,z) = = v x (§,z) . (37) 

Eq. 137(1 and a lateral average of the Stokes equa- 
tion across the computational cell i.e. (v(z)) — 

A/2 

(2/A) J v(x, z) dx, leads to {v) zz =0. As in the trans- 
verse case, the averaged velocity profile, (v), is therefore 
a linear function of z. Geometric similarity determines 
the equation for the effective slip length, namely 

L, (v(z = 0)) 

d (v(z = d)) - (v(z = 0)) ' [ ' 

Figure Ufa) represents numerical results for the normal- 
ized effective slip length, L s /d, as a function of X/d. Over 
the range shown, L s /d increases monotonically with X/d. 
As with the transverse geometry, there is no significant 




Dimensionless period X/d 




Dimensionless stripe width a/X 



FIG. 4: (Color online) (a) Normalized slip length L s /d ver- 
sus normalized pattern period X/d derived from the Stokes 
solutions for the longitudinal flow orientation. The parame- 
ter values shown are b/d = 0, 0.048 and 0.098 and a/A = 0.5. 
The straight line superimposed on the data for b/d — cor- 
responds to Eq. 11391 where L s /d~ X/d. (b) Normalized slip 
length L s /d versus normalized stripe width a/X for X/d — 1 
and o/d = and 0.098. The straight line superimposed on the 
numerical solutions corresponds to Eq. 14111 . The data points 
for b/d = 0.098 are scaled by a factor 0.5 for convenience, (c- 
d) Velocity contours corresponding to the longitudinal Stokes 
flow solutions for (c) A/d = 0.35 and (d) 10, where 6/d = 0.048 
and a/X = 0.5. The domain size is A/2 x d; the thin verti- 
cal lines designate the boundary between surfaces of no-shear 
(left) and finite slip (right). 



increase in slip length beyond X/d ~ O(10). The ab- 
solute values of L s /d are larger than in the transverse 
case. This is due to the fact that for unidirectional flow, 
the liquid above the region of no-shear always remains in 
line with the frictionless stripes and is never subject to 
any deceleration caused by the regions of finite slip. The 
functional dependence of L s /d on X/d, however, is iden- 
tical to the transverse orientation. As A/d— > 0, the slip 
length should be independent of the cell depth, d, and in- 
dependent of the type of flow (whether pressure- or shear- 
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driven). Using the analytical solutions of Philip [34 |3a| 
for the longitudinal configuration, Lauga and Stone [33| 
extracted the effective slip length for pressure driven flow 
in a cylindrical tube of radius R in the presence of alter- 
nating stripes of no-shear and no-slip (6 = 0): 



lim — = In ( 

b/R=0 R ttR \ COS 
\/R->0 X 



' X-o 



(39) 



The solutions to Eq. I|39|) , obtained by replacing the cap- 
illary radius, R, with the planar cell depth, d, are almost 
indistinguishable from the results of the full numerical 
solutions in Fig. Ufa). In this limit, the slip length L s /d 
is exactly twice that of the transverse configuration [see 
Eq. H31|l ] and scales linearly with X/d. 

An analytic expression for the effective slip, L s /d, can 
be derived in the limit X/d—>oo, by examining the flow 
field above the patterned substrates. The velocity profile 
above the no-shear surface (1) is plug-like and given by 
vi(z) — U. Above the surface of finite slip, V2(z) = 
U(z + b)/ (2b + d) . The latter result is obtained by noting 
that the shear rate, u z , is constant throughout the gap 
depth and equal to U/ (2b + d). Calculating the average 
flow speed, (u), at the upper and lower boundaries and 
substituting these into Eq. 1)38(1 leads to the expression: 



lim 

X/d^c 



d 



2X 



a b 
~ d 



(40) 



The horizontal dashed line for b/d = 0.098 and X/d > 10 
corresponds to Eq. (|40|) . 

In Fig.^Jb) is plotted the numerical solutions for L s /d 
versus a/A for X/d = 1 and b/d = and 0.098. The 
data points for b/d = 0.098 have been scaled by 0.5 for 
convenience. For b/d = and small values a/ A, a Taylor 
expansion of Eq. H39fl gives 



lim — - 

b/d=0 d 
X/d^O 



X 

ird 



(41) 



The straight line superimposed on the data in Fig. Efb) 
represents the asymptotic values given by Eq. (|41|l . The 
agreement with the analytical limit for a/A < 0.3 is very 
good. 

Velocity contours, corresponding to constant values of 
v, are shown in Figs. EJc,d). The left and right panels 
represent solutions for (c) A/d = 0.35 and (d) X/d =10 
where o/d = 0.048 and a/ A = 0.5. The vertical line de- 
notes the position corresponding to the change in bound- 
ary condition at the lower wall from no-shear (left) to 
finite slip (right). For A/d<0.35, the velocity contours 
are horizontal throughout almost the entire cell and the 
deviations from pure shear flow over a homogeneous sur- 
face are confined to a small distance from the patterned 
wall. For X/d= 10, the perturbation extends vertically 
across the cell. For large X/d, the velocity distribution 
varies from plug-like above the region of perfect slip to 
Couette-like above the region of finite slip, as assumed 
in the derivation leading to Eq. (|40|) for X/d — > oo. 



III. MD SIMULATIONS AND PARAMETER 
VALUES 

We have previously used MD simulations to investigate 
what equilibrium parameters control the degree of slip in 
simple and polymeric fluids and how the slip length de- 
pends on shear rate psl l29l|. In these previous studies, 
the wall-fluid potential was spatially homogeneous. In 
this current work, we examine the behavior of the effec- 
tive slip length for a fluid subject to planar shear in the 
presence of a heterogeneous bottom wall for the two flow 
configurations shown in Fig. ^ The wall-fluid interac- 
tions are adjusted to mimic alternating stripes of finite 
slip and no-shear by adjusting the attractive part of the 
potential to simulate more attractive and less attractive 
regions. The MD simulations described next were con- 
ducted with the LAMMPS numerical code [43j]. In what 
follows, we refer to the more attractive surface as wetting 
and the less attractive surface as non-wetting. 

The simulation cell consisted of 30720 fluid molecules 
interacting through a Lennard-Jones (LJ) potential, 



V L j(r) = 4e 



(42) 



where e and a represent the energy and length scales 
characteristic of the fluid phase. The cut-off radius was 
set to r c = 2.5 a. The parameter 6, which controls the at- 
tractive part of the potential for fluid-fluid interactions, 
was held fixed at 5 = 1. The wall-fluid (wf) parameters 
were chosen to be cr w f = 0.75 a and e w f/e = 0.8, 0.9 or 
1.0. Surfaces of finite slip in the hydrodynamic analy- 
sis corresponded to the parameter value S w { = 1.0 (i.e. 
wetting); surfaces of no-shear (or likewise perfect slip) 
corresponded to the value <5 w f = 0.1 (i.e. non- wetting) . 
For the MD simulations, we restricted our study to the 
case a/A = 1/2 such that the wetting and non- wetting 
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FIG. 5: (Color online) Average normalized fluid density, 
p(z)a 3 , above the wetting (<5 w f = 1-0: triangles) and non- 
wetting stripes (5 w f = 0.1: circles). The parameter values 
shown are (a) a = 1.04<r and (b) o= 133.48 a for e wf /e = 0.8. 
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FIG. 6: (Color online) Average normalized velocity profile, 
(u) r/a corresponding to the transverse flow configuration for 
£ wf /£ = 0.8 and a/ a = 1.04, 4.17 and 133.48. The nearly hori- 
zontal profile shown in (4), which designates a plug- like profile 
representative of surfaces with large slip lengths, was obtained 
by setting the wall-fluid potential parameter to 5 w f =0.1 
along both walls. This choice effectively reproduces non- 
wetting bounding walls. 



portions of the substrate occupy equal areas. 

The upper and lower walls of the simulation cell each 
consisted of 12288 molecules distributed between two 
(111) planes of an FCC lattice with density p w = 4p, 
where p = 0.81 <r~ 3 is the density of the fluid phase. The 
fluid was confined to a fixed height d = 20.15 tr; the cell 
volume was 266.96 a x 7.22 a x d for the transverse ge- 
ometry. To eliminate any finite size effects for the lon- 
gitudinal geometry, the system size along the y-axis was 
doubled in length to 14.45 er, requiring simulations with 
61440 fluid molecules. For either configuration, periodic 
BCs were enforced along the x and y axes. The fluid was 
held at a constant temperature T = l.le/fcs by means 
of a Langevin thermostat 0] with a friction coefficient 
t _1 . Here, ks is the Boltzmann constant. This damping 
term is only applied to the coordinate equation perpen- 
dicular to the direction of flow |28l l36| . The equations 
of motion were integrated using the Verlet algorithm 
with a time step At — 0.005 t, where r = y/ ma 2 /e rep- 
resents the characteristic time set by the LJ potential 
and m is the monomer mass. The fluid was subject to 
steady planar shear by translating the upper wall at a 
constant speed U; the lower, patterned wall remained 
stationary. In all the simulations, the speed of the upper 
wall was held fixed at [/ = 0.5cr/r. After an equilibration 
period exceeding 10 r, the fluid velocity profile was ob- 
tained by averaging the instantaneous monomer speeds 
in slices Az = 0.1 a for a time interval At « 3 • 10 4 r. 
The Reynolds number, based on the upper wall speed U, 
the wall separation d and the fluid shear viscosity (de- 
termined previously [HH^ to be = 2.2 ± 0.2er/cr 3 
for comparable shear rates) was estimated to range from 
2 — 5, indicative of negligible inertial effects and laminar 
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FIG. 7: (Color online) Average normalized velocity profile, 
{v)t/ct, corresponding to the longitudinal flow configuration 
for £ wf /e = 0.8 and a/ a = 1.04, 4.17 and 133.48. 

flow conditions. In fact, this estimate provides only an 
upper bound on the Reynolds number, since the actual 
fluid velocity for surfaces comprising regions of finite and 
infinite slip is significantly smaller than the upper wall 
speed. In our studies, use of the fluid flow speed further 
reduces Re by a factor of up to 2. We conclude that 
the small Reynolds numbers characterizing the MD sim- 
ulations are consistent with the theoretical restriction for 
the Stokes flow solutions obtained in the limit Re = 0. We 
also note that the numerical solutions to the Stokes equa- 
tion for the longitudinal geometry are valid irrespective 
of the value of the Reynolds number because the unidirec- 
tional flow causes the inertial term in the Navier-Stokcs 
equation to vanish identically. 

IV. RESULTS OF MD SIMULATIONS FOR 
TRANSVERSE AND LONGITUDINAL FLOW 

The two sets of curves in Fig. 03 show the average 
normalized fluid density, p(z)a 3 , for the transverse flow 
configuration in the region above the wetting and non- 
wetting stripes for e w f/e = 0.8 and a/a — 1.04 and 133.48. 
The choice a — 133.48 a represents the accommodation of 
only two stripes at the lower wall within the Couette cell. 
The oscillations near the upper and lower boundaries re- 
flect the molecular layering caused by the presence of 
dense walls [3(j. Increasing the attractive part of the 
LJ potential generates larger peak maxima and more os- 
cillations. Above either type surface, the density oscilla- 
tions persist for about 4 — 6 molecular diameters from the 
wall. Decreasing the strength of the attractive interac- 
tion shifts the first peak maximum away from the lower 
wall. Also, the fluid density above the wetting stripes 
is found to increase with a/ a. The density profiles cor- 
responding to longitudinal flow configurations are quite 
similar to the ones shown here. 

Figure [S] shows representative velocity profiles across 
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FIG. 8: (Color online) Comparison of the effective slip length, 
L s , as extracted from the MD simulations (symbols) , with nu- 
merical solutions of the Stokes equation (continuous lines) for 
transverse flow. The local slip length, b, as extracted from the 
MD simulations, decreases with increasing wall-fluid attrac- 
tion energy, namely b/a = 1.97, 1.36 and 0.95 for e w f/e = 0.8 
(o), 0.9 (o) and 1.0 (V), respectively. The dashed horizontal 
lines for a/a > 100 correspond to Eq. I3UI . Inset: MD re- 
sults showing collapse of the effective slip length L 3 /a when 
rescaled by the quantity e^ f , versus a/a. 



the cell depth for transverse flow with e w i/e — 0.8 and 
a/a = 1.04, 4.17 and 133.48. Shown for comparison 
is the velocity profile corresponding to the case of uni- 
formly no n- wetting walls where 5 w f = 0.1 holds for both 
surfaces. Decreasing the wall-fluid interaction leads to 
a high degree of slip and a plug-like velocity field. The 
remaining three profiles increase linearly with z/ a, as ex- 
pected for a fluid subject to planar shear, except in the 
vicinity of the lower wall. Significant deviations from 
linearity occur for large stripe widths. These oscillations 
are caused by the difference in the positions of the fluid 
density maxima above the wetting and non-wetting re- 
gions [see Fig.|3Jb)]. As evident from the velocity profile, 
the degree of slip increases with increasing values of a. 

Figure 13 shows the computed velocity profiles for lon- 
gitudinal configurations. The behavior is similar to that 
shown in Fig. for the transverse orientation, but the 
amplitude of the oscillations near the lower wall is sig- 
nificantly larger. In this case, the degree of slip does not 
increase monotonically with a. The smallest stripe width 
generates the second largest slip velocity in Fig. [7| As 
the stripe width increases, it is found that the wetting re- 
gions induce stronger molecular ordering in the first fluid 
layer adjacent to the wall, causing a reduction in the slip 
length, as noted in Fig. ^3 

For direct comparison to the hydrodynamic predic- 
tions, it was necessary to extract the actual values of the 
local slip length, b, representative of the surfaces char- 
acterized by <5 w f = 1.0, for input values to the boundary 
conditions used in computing the solutions to the Stokes 
equation. This was accomplished in the MD simulations 



FIG. 9: (Color online) Direct comparison of the effective slip 
length extracted from the MD simulations (symbols) with the 
numerical solutions of the Stokes equation (continuous lines) 
for longitudinal flow. The local slip length, b, as extracted 
from the MD simulations, varies with the LJ wall-fluid inter- 
action energy, e w f as 6/a = 1.97, 1.36 and 0.95 for e w f/e = 0.8 
(o), 0.9 (o) and 1.0 (V), respectively. The local slip lengths 
are observed to be independent of the flow orientation. 



by extrapolating the average velocity profile at the top 
wall to a speed U for different values of a imposed on 
the lower wall. The extrapolated distance b was found to 
depend on the wall-fluid interaction energy but not the 
shear rate in the fluid nor the flow orientation. As ex- 
pected, the values of b decreased with increasing value of 
the wall-fluid interaction energy, namely b/a— 1.97, 1.36 
and 0.95 for e w f/e = 0.8, 0.9 and 1.0, respectively. By 
contrast, the local slip length for the flat velocity profile 
shown in Fig.Elfor uniformly non- wetting walls was found 
to be (362±10) a. Given that this slip length significantly 
exceeds the wall separation, the choice (5 w f = 0.1 approx- 
imates very well the behavior of surfaces of perfect slip 
(i.e. no-shear) assumed in the continuum calculations. 

The composite or effective slip length, L s , was deter- 
mined in the MD simulations by linear extrapolation be- 
low the stationary lower surface of the velocity profile to 
the value zero. Figure |S] represents a plot of L s /a with 
increasing normalized stripe width, a/a, and increasing 
wall-fluid interaction strength, e w f, for transverse flow 
configurations. The MD results (symbols) show a sharp 
increase in slip length for a/ a < 10 and saturation to a 
constant value beyond a/a > 100. 



V. DISCUSSION 

As described in Section III Al and for fixed values of 
a/ A, the effective slip length derived from hydrodynamic 
considerations depends only on the ratios X/d and b/d. 
The molecular length scale, a, plays no part in the anal- 
ysis. For direct comparison to the MD results, it was 
therefore necessary to multiply the numerical values of 
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L s /d, X/d and b/d from the Stokes solutions with the 
value of the wall separation, d = 20.15 a, used in the 
MD simulations. The largest ratio, a/d = 6.62, acces- 
sible to the MD simulations was only limited by com- 
putational resources. The solid lines shown in Fig. |S| 
represent solutions of the Stokes flow equation for trans- 
verse flow, as discussed in Section Til Al The agreement 
between the continuum predictions and the MD simula- 
tions is excellent for a /a > O(10); significant deviations 
occur for a/a < 0(1). The asymptotic predictions given 
by Eq. l|5U)l for X/d = 2a /d — > oo are designated by the 
dashed horizontal lines in Fig. [S] 

The Green-Kubo type analysis of Barrat and Boc- 
quet [23l for homogeneous surfaces characterized by 
a single wall-fluid interaction energy predicts that the 
slip length scales as e~f provided the in-plane structure 
factor, fluid contact density and in-plane diffusion coef- 
ficient characteristic of the first fluid layer remain rela- 
tively constant. The results shown in the inset of Fig. El 
for the transverse geometry confirm this prediction for 
the range a < 10 a, even for the case of a composite 
potential where the wall-fluid interaction alternates be- 
tween two values of S w i . This collapse fails above a > 10 a 
where the continuum solutions show excellent agreement 
with the molecular simulations. This behavior suggests 
that for a/a < O(10), the effective slip length is mostly 
determined by the molecular scale frictional properties 
between the first fluid layer and the lower wall. For 
a / a <L O(10), however, the effective slip length is set 
by the wall separation d, the pattern lengthscales a and 
A and the local slip length b. The transition region 
8 0,/a < 30 therefore contains mesoscopic information 
from both the molecular and hydrodynamic descriptions. 

The deviation between the MD simulations and the 
Stokes solutions below a/a < O(10) can be understood 
as follows. The lower wall is comprised of a potential 
whose interaction strength alternates between wetting 
and non-wetting values with a periodicity set by the 
stripe width a, which approaches the molecular scale. 
The fluid molecules no longer experience uninterrupted 
stretches of wetting and non- wetting regions; instead, the 
fluid molecules are exposed to an effectively roughened 
surface with molecular scale corrugations. These corru- 
gations trap the fluid molecules, thereby suppressing slip 
at the wall-fluid interface. The commensurability be- 
tween the fluid molecular size and the wall corrugation 
size can in fact lead to a no-slip condition for slightly 
larger values of e w f 01 • It is therefore not surprising that 
the effective slip length for the transverse configuration, 
as shown in Fig. EI decreases sharply with decreasing val- 
ues of a. This effect also explains why for the smallest 
values of a the slip length L s is even smaller smaller than 
the local slip length obtained for a fluid confined between 
two identical walls both characterized by the same value 
(5 w f = 1.0. For example, for e w f/e = 0.8 and a/a = 1.04, 
we find that b/ a = 1.97 but L s /a = 0.h\ 

Figure shows the behavior of L s /a as a function of 
stripe width, a/ a, and increasing interaction strength, 



e w f, for longitudinal flow configurations. The results of 
the MD simulations (symbols) show a sharp decrease in 
slip length below a/a ^ 10 followed by a rapid rise. The 
effective slip lengths have similar magnitudes for very 
small and very large values of a. Once again, there is ex- 
cellent agreement between the Stokes flow solutions and 
the MD simulations for a/a > O(10) but strong devi- 
ations below this value. In contrast to the transverse 
configuration, however, the MD results predict much 
larger effective slips than the continuum solutions for 
a/a < O(10). Because of the translational invariance of 
the flow inherent in this case, the molecular scale rough- 
ness set by the composite potential at the bottom wall 
cannot diminish the slip length. The reduction in molec- 
ular ordering above the wetting regions with decreasing 
stripe width, however, leads to an increase in the slip 
length which exceeds the slip lengths obtained for the 
transverse configuration as well as the continuum predic- 
tions. 

Previous MD simulations of Newtonian and non- 
Newtonian fluids have demonstrated that the slip length 
for surfaces characterized by a single wall-fluid poten- 
tial correlates strongly with the degree of molecular 
ordering in the first fluid layer adjacent to the wall 
@ mil |H HH H3. The more orderly the molec- 
ular organization, as reflected by the maximum value of 
the in-plane structure function, S m&Xl the smaller the slip 
length. To test these predictions for the case of patterned 
walls in a longitudinal orientation, we computed the max- 
imum value of the in-plane structure function within the 
first fluid layer above the wetting and non-wetting re- 




FIG. 10: (Color online) The dominant peak in the in-plane 
fluid structure factor evaluated separately above the wetting 
(o) and non-wetting (o) regions for longitudinal flow and 
£wf/e = 0.8. Inset: MD results showing the strong correla- 
tion between the effective slip length L 3 /a (circles: data from 
Fig. |5J and the quantity A(S' max p c (j 2 ) _1 (triangles), which 
characterizes the degree of molecular ordering within the first 
fluid layer above the wetting stripes. The value of the fitting 
parameter A is 92.4. 
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was estimated from the position of the first minimum in 
the density profile above a wetting stripe. The contact 
density p c was identified with the maximum of the fluid 
density within the first fluid layer. The structure func- 
tion was computed according to S(q)= | Y^i* e iqv \ 2 /Ne, 
where Ne, is the number of molecules in the first fluid layer 
adjacent to either a wetting or non-wetting surface. As 
shown in Fig. ^3 the molecular ordering adjacent to a 
wetting region is far stronger and increases with increas- 
ing stripe width, a. By contrast, the molecular ordering 
adjacent to the non-wetting region is unaffected by the 
stripe width, a, except for the smallest value shown. The 
inset in Fig. ^J] demonstrates the correlation between the 
effective length and the parameter (SmaxPcf 2 ) -1 as es- 
timated above the wetting regions. Here, the values of 
L s /<j versus a/a from Fig. are plotted alongside the 
quantity A(S max p c a 2 )^ 1 , where A — 92A is a fitting pa- 
rameter. In the limit a/a<O(l0), the strong correlation 
between L s and A(S max p c a 2 )~ 1 establishes that the in- 
crease in effective slip length for narrow stripe widths 
is mainly caused by the reduction in molecular ordering 
within the first fluid layer above the wetting zones. 

The BCs used in the continuum analysis correspond to 
stripes of finite (or no) slip and no shear (i.e. b = oo). 
We repeated the analysis in Section II by replacing the 
no-shear BC with a second slip BC to define surfaces 
comprising alternating stripes of small (b/a — 1.97) and 
large slip (b/a = 362, as extracted from case (4) shown 
in Fig. 6). For the transverse configuration, the curve 
corresponding to b/a — 1.97 in Fig. 8 showed a slight 
decrease in L s of about 3% for a/a > 30, whereas the 
longitudinal configuration generated a decrease of up to 
9% with respect to the values shown in Fig. |5] 

VI. SUMMARY 

We have investigated the behavior of the slip length 
in Newtonian liquids subject to planar shear in a Cou- 
ette cell with mixed surface boundary conditions. The 
upper wall is modelled as a homogenous surface with fi- 
nite or no-slip moving at a constant speed above a lower 
stationary wall patterned with alternating stripes repre- 
senting regions of no-shear and finite or no-slip. The ve- 
locity fields and effective slip lengths are computed both 
from molecular dynamics (MD) simulations and solution 
of the Stokes equation for flow parallel (longitudinal case) 
or perpendicular (transverse case) to the stripe pattern. 
Detailed comparison between the results of the hydro- 
dynamic calculations and MD simulations shows excel- 
lent agreement when the length scale of the substrate 
pattern geometry is larger than O(10 a), where a de- 
notes the fluid molecular diameter as set by the Lennard- 



Jones interaction. The effective slip length then increases 
monotonically with a/ a to a saturation value. For the 
transverse case, the Stokes flow solutions predict an effec- 
tive slip larger than the MD results when a/a ~ O(10). 
This discrepancy is understood from a molecular point 
of view since a narrowing of the regions subject either 
to no-shear or no-slip essentially establishes a roughened 
surface. The molecular scale corrugation created by the 
composite wall potential strongly reduces the effective 
slip length below the hydrodynamic results. This surface 
roughening effect is not present for the longitudinal flow 
configuration since the fluid molecules are transported 
along homogeneous stripes representing regions of cither 
no-shear or finite slip. In this case, however, the 2D 
fluid structure factor above the non-wetting stripes (re- 
gions of perfect slip or equivalently no-shear) decreases 
for a/ a < 0(10), which enhances the effective slip lengths 
above the values predicted by the hydrodynamic solu- 
tions. On the molecular level, the strong correlation ob- 
served between the effective slip length and the product 
(pc&nax) confirms that a reduction in molecular order- 
ing within the first fluid layer generates an increase in the 
effective slip length. 

Detailed comparison between continuum computations 
and molecular dynamics simulations is of increasing im- 
portance to the development of hybrid computational 
schemes 0, 0, lE^- These algorithms are de- 
signed to stitch together hydrodynamic solutions ob- 
tained from continuum equations with the molecular 
scale solutions obtained from MD simulations or other 
microscopic solvers. It has been demonstrated that the 
spatial coupling across this wide range in length scales 
can be achieved by implementation of constraint dynam- 
ics within an overlap region. We hope that our studies 
of shear driven flow along surfaces with mixed boundary 
conditions will complement ongoing efforts using hybrid 
codes. The system and results described here offer an 
interesting test case for better understanding of the in- 
termediate region bridging the behavior of fluids from the 
nanoscale to microscale dimensions. 
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